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In this paper, we investigate numerically the flow of an electrically conducting fluid in a cylindrical 
Taylor- Couette flow when an axial magnetic field is applied. To minimize Ekman recirculation due 
to vertical no-slip boundaries, two independently rotating rings are used at the vertical endcaps. 
This configuration reproduces setup used in laboratory experiments aiming to observe the Magne- 
toRotational Instability (MRI). Our 3D global simulations show that the nature of the bifurcation, 
the non-linear saturation, and the structure of axisymmetric MRI modes are significantly affected 
by the presence of boundaries. In addition, large scale non-axisymmetric modes are obtained when 
the applied field is sufficiently strong. We show that these modes are related to Kelvin-Helmoltz 
destabilization of a free Shercliff shear layer created by the combined action of the applied field 
and the rotating rings at the endcaps. Finally, we compare our numerical simulations to recent 
experimental results obtained in the Princeton MRI experiment. 
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I. INTRODUCTION 

The magnetorotational instability has provided a 
simple explanation of the longstanding problem of 
rapid angular momentum transport in accretion disks 
around stars and black holes Balbus and Hawley, 
rediscovering an instability first studied by Velikhov 
Q and Chandrasekhar have shown that keplerian 
flows of accretion discs, for which the Rayleigh criterion 
predicts axisymmetric hydrodynamical stability, can be 
destabilized in the presence of a magnetic field. More 
precisely, the MRI is a linear instability occurring when 
a weak magnetic field is applied to a rotating electri- 
cally conducting fluid for which the angular velocity 
decreases with the distance from the rotation axis, i.e. 
9r(S1 2 ) < 0. Linear stability analysis indicates that the 
most unstable mode is axisymmetric and associated with 
a strong radial outflow of angular momentum. Non- 
linear evolution of the MRI is of primary importance, 
since saturation of the instability eventually yields a 
magnetohydrodynamical turbulent state, enhancing the 
angular momentum transport. 

During the last decade, there has been a lot of effort 
to observe the MRI in a laboratory experiment. To this 
end, most groups use a magnetized Taylor- Couette flow, 
i.e. the viscous flow of an electrically conducting fluid 
confined between two differentially rotating concentric 
cylinders, in the presence of an externally imposed mag- 
netic field. For infinitely long cylinders, the ideal laminar 
Couette solution is given by: 

m=M + ^ a) 

in which A\ = (VL2t\ — ^l\r\)j{r\ — r\) and A 2 = 
r\r — ^2)/ {r\ — rf), Qi and ^2 are respectively the 
angular velocity of inner and outer cylinder, and 7*1, r 2 



are the corresponding radii. The Rayleigh criterion pre- 
dicts axisymmetric linear stability if ^2/^1 > (V1/V2) 2 , 
ensuring that the specific angular momentum increases 
outward. However, MRI may arise if ^2/^1 < 1 and 
non-ideal MHD effects are small Q. 

Several experiments are currently working on MRI. 
The Princeton experiment has been designed to ob- 
serve this instability in a Taylor- Couette flow of liquid 
Gallium, with an axial applied magnetic field |5j. So 
far, MRI has not been identified, but non-axisymmetric 
modes have been observed when a strong magnetic field 
is imposed [6|. The PROMISE experiment, in Dresden, 
is based on a similar set-up, except that the applied field 
possesses an azimuthal component. Axisymmetric trav- 
eling waves have been obtained, and identified as being 
Helical MRI, an inductionless instability different from 
but connected to the standard MRI 0. A few years ago, 
MRI was claimed to have been obtained in a spherical 
Couette flow of liquid sodium (hereafter "Maryland 
experiment") [8]. In this experiment, in which the outer 
sphere was at rest and a poloidal magnetic field was 
applied, non-axisymmetric oscillations of both velocity 
and magnetic fields were observed, together with an 
increase of the torque on the inner sphere. However, 
it has been shown recently that these oscillations were 
more likely related to instabilities of magnetic free shear 
layers in the flow @, PlJ- 

Experimental observation of the MRI is considerably 
complicated by the presence of vertical boundaries. In- 
deed, no-slip boundary conditions at the vertical end- 
caps (for instance rigidly rotating with either the outer 
or the inner cylinder) induce an imbalance between the 
pressure gradient and the centrifugal force, and drive a 
meridional Ekman recirculation in addition to the az- 
imuthal Couette flow. An inward boundary-layer flow 
near the endcaps is balanced by strong outward jet at the 



2 



midplane. This Ekman recirculation has two important 
consequences for laboratory MRI: first, the meridional 
flow transports angular momentum outward, decreasing 
the free energy available to excite MRI. Moreover, the 
Ekman flow is a source of hydrodynamical fluctuations, 
beginning at Reynolds numbers > 400 with oscillations 
of the radial jet PJJ, which strongly complicates the iden- 
tification of MRI modes. 

In the Princeton MRI experiment, this problem has 
been circumvented by replacing the rigid endcaps at 
the top and the bottom by two rings that are driven 
independently, reducing the imbalance due to viscous 
stress. It has been shown that a quasi- keplerian flow 
profile (ri/r2) 2 < ^2/^1 < 1 in a short Taylor- Couette 
cell can be kept effectively stable up to Re ~ 10 6 if 
appropriate ring speeds are used [5]. Similarly, the 
PROMISE group found that split rings (attached to 
the cylinders rather than independently driven) led to a 
significant reduction of the Ekman pumping and a much 
clearer identification of the HMRI [12j . 

In this article, we report three-dimensional numeri- 
cal simulations inspired by these experimental configura- 
tions, especially the Princeton MRI experiment. In the 
first section, we present the equations and the numerical 
method used. In section II, we study how the structure 
and the saturation of the magnetorotational instability 
arising in Taylor- Couette flow is modified by the presence 
of vertical boundaries. In section III, we report new non- 
axisymmetric instabilities generated by a free shear layer 
in the flow when the applied field is sufficiently strong 
compared to the rotation. Finally, a comparison with 
results from the Princeton MRI experiment is presented. 



II. EQUATIONS 

We consider the flow of an electrically conducting fluid 
between two co-rotating finite cylinders. r\ is the radius 
of the inner cylinder, r2 = 3ri is the radius of the outer 
cylinder, and d = T2 — v\ is the gap between cylinders, 
f^i and 0,2 are respectively the angular speed of the inner 
and the outer cylinder. The height of the cylinders H is 
fixed such that we have the aspect ratio Y = H/d = 2, as 
in the Princeton MRI experiment. A uniform background 
magnetic field Bo = B§e z is imposed by external coils. 

The governing equations for this problem are the 
Navier- Stokes equations coupled to the induction equa- 
tion : 



p— + p(u-V)u 



-VP + pV 2 w + j X B . (2) 



f) 1 

= V X (« X B) + —V 2 B . (3) 

at pocr 

where p is the density, v the kinematic viscosity, 
r] = 1/(07x0) is the magnetic diffusivity, u is the fluid 



velocity, B the magnetic field, and j = /Xq x V X B 
the electrical current density. The problem is charac- 
terized by three dimensionless numbers: the Reynolds 
number Re = (VLirid)/v, the magnetic Reynolds 
number Rm = (Clirid)/rj and the Lundquist num- 
ber S = Bori/(y/p/jLov)' Alternatively, the applied 
field can be measured by the Hartmann num- 
ber M = BqVi/ y/p/jLorjv, or the Elsasser number 
A = Bq/(ppqT]Q). It is also useful to define the magnetic 
Prandtl number Pm = Z//77, which is simply the ratio 
between Reynolds numbers. At the top and bottom, the 
endcaps are divided at r = (7*1 + 7*2) /2 between two in- 
dependently concentric rotating rings. Rotation rates of 
inner and outer rings are given respectively by £7 3 and ^4. 

These equations are integrated with the HERACLES 
code [13]. Originally developed for radiative astrophys- 
ical compressible and ideal-MHD flows, HERACLES 
relies on a finite volume Godunov method. The code has 
been parallelized with the MPI library and implemented 
in Cartesian, cylindrical and spherical coordinates. 
We have modified this code to include fluid viscosity, 
magnetic resistivity, and suitable boundary conditions 
for velocity and magnetic fields. Note that HERACLES 
is a compressible code, whereas laboratory experiments 
generally involve almost incompressible liquid metals. In 
fact, incompressibility corresponds to an idealization in 
the limit of infinitely small Mach number (Ma). In the 
simulations reported here, we used an isothermal equa- 
tion of state with a sound speed such that Ma < 0.1, 
following the approach of [H[,[l5|. Typical resolutions 
used in the simulations reported in this article are 
(Nr.N^.Nz) = [200,64,400]. For some runs, we have 
checked that doubling this resolution does not affect 
results. 

Except where indicated, the Reynolds number has 
been fixed at Re = 4000, whereas Rm and S vary widely. 
On the cylinder sidewalls, at r = r\ and r = 7*2, we im- 
pose perfectly electrically conducting boundaries : 



B r = , dB z /dr = , dB^/dr = (4) 

At the vertical endcaps, we use the so-called pseudo- 
vacuum boundary conditions, namely the magnetic field 
is forced to be normal to the boundary: 



B r = 



B 6 = 



dB z /dz = 



(5) 



With these conditions, no magnetic torques are exerted 
through the boundaries. For the velocity field, no-slip 
boundary conditions are used for radial and axial com- 
ponents. The angular velocity matches the rotation rates 
of cylinders and rotating rings at radial and axial bound- 
aries. 
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III. MAGNETO-ROTATIONAL INSTABILITY 
IN FINITE GEOMETRY 

A. Purely hydrodynamical simulations 

Accretion discs result from the balance between gravi- 
tation and centrifugal force, leading to an angular veloc- 
ity profile such that ft ~ r -3 / 2 . Since the Rayleigh cri- 
terion states that any inviscid flow with rotation profile 
ft ~ r 6 is linearly stable to the centrifugal instability if 
S > — 2, accretion discs are expected to be hydrodynami- 
cally linearly stable. Although Keplerian profiles can not 
be reproduced in the laboratory, it is possible to drive 
a Taylor- Couette flow in the so-called quasi-keplerian 
regime, i.e. with an angular momentum which increases 
outward to satisfy Rayleigh stability criterion, but which 
is MRI unstable (the angular velocity decreases outward). 
Except where indicated, all simulations reported here use 
the following rotation rates for the cylinders and rings: 

«i = 1. , ft 2 = 0.13 , ft 3 = 0.45 , Ol = 0.16 (6) 

where rotation rates have been normalized by fl±. The 
ratio f^/f^i ensures a quasi-keplerian regime, whereas 
and have been carefully chosen to strongly reduce 
Ekman recirculation and therefore reinforce the hydro- 
dynamical stability of the flow [l5| . 

Figure [1] shows the flow obtained with the aforemen- 
tioned standard parameters but without magnetic field, 
and for an axisymmetric simulation. Note that the ra- 
dial profile of angular velocity matches almost perfectly 
the theoretical Couette solution (Fig. [TJ top-left) that 
would be obtained with infinite or ^-periodic cylinders. 
In fact, figure [TJ top-right shows that the Couette profile 
prevails not only in the midplane, but in most of the do- 
main. As expected, the use of rotating rings at the end- 
caps strongly reduces the Ekman flow, dividing the usual 
two big Ekman cells into eight smaller and less intense 
cells (Fig. [IJbottom- left). This structure is unsteady and 
the corresponding weak poloidal flow fluctuates as time 
evolves. The amplitudes of u r and u z are one order of 
magnitude smaller than for endcaps rotating with the 
outer cylinder. 

To compare our numerical results to theoretical pre- 
dictions, it is useful to define a local shear parameter q: 
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(7) 



For instance, in a keplerian flow q would have the 
constant value —3/2. Figure [TJbottom-right shows 
q in the (r, z)-plane for our system. In the bulk of 
the flow, q is between —1.5 and —2 and matches the 
value corresponding to an ideal circular Couette flow. 
However, close to the endcaps, a very strong shear 
is generated where the rings meets, characterized by 
q < —3.5, reflecting Rayleigh-unstable flow. This 
corresponds to a cylindrical Stewartson free shear layer 
[l6| at r = (ri + r2)/2, created by the jump of the 
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FIG. 1: Purely hydrodynamical axisymmetric 
Taylor- Couette flow, obtained for Re = 4000, Rm = 15, 

S = and rotations rates given by eq.flBJ). Top, left: 
profile of ft at z = h/4. Top, right: angular velocity ft 

in (r, z) plane. Bottom, left: U z in the (r, z) plane 
(streamlines indicate meridional flow). Bottom, right: 
Shear q (see eq. ([7]) ) in the (r, z) plane. The use of 
rings at endcaps allows the flow to be very close to 
theoretical circular Couette flow. There is a weak 
poloidal recirculation, with 8 unsteady cells. Except 
close to the boundaries, the shear is relatively uniform. 



angular velocity between the two independently rotating 
rings. See [17j for a more detailed numerical study of 
this layer. Although one could expect a destabilization 
of this shear layer to large scale non- axisymmetric modes 
[l8| , this has not been observed in our simulations. Note 
that the layer does not extend very far from the endcaps, 
suggesting that it is disrupted by small-scale instabilities. 



B. Imperfect bifurcation of the magnetorotational 
instability 

Let us now consider the magnetized problem, meaning 
that a homogeneous axial magnetic field is now applied 
at the beginning of the simulation. As expected from 
global and local linear analyses, the flow is destabilized 
to MRI for sufficiently large Rm and S. Figure [2] shows 
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FIG. 2: Stability boundaries in the plane of magnetic 
Reynolds number, Rm, versus Lundquist number, S. 
Black curve indicates the marginal stability curve of the 
MagnetoRotational Instability (MRI) obtained by 
interpolation from 2D simulations. Red curve is the 
marginal stability curve of non-axisymmetric modes 
generated by Kelvin-Helmoltz destabilization of the 

MHD detached shear layer, obtained from 3D 
simulations (black dots). Dashed square indicates the 
parameter space accessible to the Princeton experiment. 



the marginal stability curve of the MRI (black curve) in 
the plane defined by these parameters. In the domain 
explored here, MRI modes are always axisymmetric. 
Linear stability analysis predicts non-axisymmetric MRI 
modes for Rm > 50, which is larger than considered 
here. Note that the non-axisymmetric modes in Fig. [2] 
are not standard MRI but are related to the instability 
described in the next section. If the applied field is too 
strong, MRI modes are suppressed, and the flow is resta- 
bilized to a laminar state dominated by the azimuthal 
velocity. As a consequence, MRI modes are unstable 
only within the pocket-shaped interior of the black curve. 

Figure [3] shows the structure of the MRI mode obtained 
in the saturated regime, for Rm = 15 and S = 2.3. 
Here, the structure of the MagnetoRotational Instabil- 
ity mainly consists of two large poloidal cells, as shown 
in the left panel. The corresponding radial magnetic 
field is shown on the right. At the endcaps, the fluid is 
strongly ejected and generates an outflowing narrow jet. 
The recirculation takes the form of a broad inflowing jet 
[l5| . Consequently, we see that the MRI modes resem- 
ble the hydrodynamic Ekman circulation obtained when 
the endcaps corotate with the outer cylinder, except that 
the circulation is reversed. In fact, MRI modes are rela- 



FIG. 3: Structure of the MagnetoRotational Instability 
(MRI) mode when a magnetic field is applied in the 
z-direction, obtained for Rm = 15 and S = 2.3. Left: 
Radial velocity U r in the (r, z) plane. Curves are 
streamlines of the poloidal flow. Right: radial magnetic 
field component i? r (normalized by the applied field B%) 
with magnetic field lines contained in the (r, z) plane. 
For these boundary rotations [eq. MRI is 
characterized by a strong inflowing jet. 



tively similar to the classical Taylor vortices obtained in 
hydrodynamical Taylor- Couette flow. 

To follow the evolution of this large scale MRI mode, 

we compute the global quantity A = \Jv~ x J v B^dV, 

where V is the total volume between cylinders. Figure 
2] shows the saturated value of A as a function of the 
Lundquist number, for Rm = 15. To understand how 
the MRI is generated in our finite geometry, it is useful 
to compare this bifurcation diagram to results obtained 
with periodic (or infinite) cylinders in the z direction. In 
Figure |4j the black curve thus indicates results obtained 
with our no-slip rotating endcaps, whereas the red curve 
shows the diagram obtained when vertically periodic 
boundary conditions are used. 

With periodic boundary conditions, the bifurcation of 
A is characterized by a well defined critical Lundquist 
number S c above which MRI is observed (S c ~ 0.25 for 
Rm = 15). The magnetic field follows a classical square 
root law A ~ ±V S — S c slightly above the MRI thresh- 
old S c , and is zero below S c . In this case, the solution is 
relatively similar to the one obtained with no-slip bound- 
aries, except that the axial position of the outflowing jet 
is made arbitrary by the periodic boundaries, and the 
linear MRI modes are sinusoidal in the z direction. In 
the linear regime, MRI modes correspond to the unsta- 
ble branch appearing from the coalescence of two com- 
plex conjugate MagnetoCoriolis (MC) waves @. If these 
MC-waves are linearly stable (as in our simulations), the 
non-linear transition to MRI corresponds, at least close 
to the threshold S c , to a classical supercritical pitchfork 
bifurcation: 



A = ii A - A 3 



(8) 
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FIG. 4: Top: Bifurcation diagram of the radial 
magnetic energy A as a function of S for Rm = 15. Red 
curve correspond to periodic boundary conditions in z, 
whereas the black one corresponds to results obtained 
with no-slip rotating rings at the endcaps. In the latter 

case, the bifurcation is an imperfect supercritical 
pitchfork bifurcation due to the poloidal flow. Bottom: 
Schematic representation of perfect (red) and imperfect 
(black) pitchfork bifurcations of the MRI. Each branch 
corresponds to a different structure of the MRI mode: 
insets show MRI mode obtained from rotation rates 
given by eq.flS) (top- inset) or eq. (fTU|) (bottom- inset). 



where the dot denotes time derivative, and [i oc S— S c is a 
control parameter. Each of the two supercritical branches 
of solutions (A = ±y/]T) is a transposition of the other 
corresponding to a phase-shift of one-half period, reflect- 
ing the arbitrary position of the MRI jet. For instance, 
the upper branch will correspond to an MRI mode with 
an outflow in the midplane, whereas the lower branch 
will represent MRI mode with inflow in the midplane. A 
schematic view of this supercritical bifurcation is shown 
in the lower panel of Fig. |4]by the red line. 

However, no-slip vertical boundary conditions signifi- 



cantly change the nature of the bifurcation: with finite 
cylinders, we see in figure SJtop (black curve) that the 
evolution of A is more gradual, and the definition of a 
critical onset S c is not well defined. In a finite geometry, 
the magnetorotational instability rather corresponds to 
an imperfect supercritical pitchfork bifurcation: 

A = fiA - A 3 + h (9) 

Here, h is a symmetry-breaking parameter reflecting the 
absence of z-periodicity of the solutions, and forbidding 
the A —> —A symmetry obtained in the case of periodic 
or infinite cylinders. Figure HJbottom shows a schematic 
diagram of supercritical pitchfork bifurcations, in perfect 
and imperfect cases. Physically, this symmetry breaking 
means that for any value of the applied magnetic field, 
there is always a non-zero solution corresponding to 
the small poloidal recirculation created by the endcaps 
(see Figure [TJbottom-left). Therefore, with realistic 
boundaries, the magnetorotational instability continu- 
ously arises from this residual magnetized return flow, 
so that the axial position of the MRI outflowing jet is no 
longer arbitrary. This can be regarded as an analogue 
of the transition encountered in finite Taylor- Couette 
flow, in which Taylor vortices gradually emerges from 
the Ekman flow located at the endcaps [19]. Figure 
[4]-top also shows that for strong applied field, the 
magnetic tension suppresses MRI modes. It is inter- 
esting to note that again, this restabilization is much 
more gradual in finite geometry than in the periodic case. 

The modified nature of the bifurcation in the presence 
of realistic boundaries has important implications for 
experimental observations of the instability, because 
it strongly complicates the definition of an onset for 
the MRI. Since a recirculation is always generated 
by the boundaries, it will be difficult to distinguish 
an MRI mode from a simple magnetized meridional 
return flow, at least close to the instability threshold. 
If the generation of the magnetorotational instability 
in laboratory experiments corresponds to a strongly 
imperfect pitchfork bifurcation, as suggested by our 
numerical results, one can expect the following: 

- First, the MRI will appear with the imperfect scaling 
obtained after integration of equation (j9j) instead of a 
classical square root law. On the one hand, this means 
that the transition to instability in bounded systems 
will not be as clear-cut as formerly expected, especially 
in experiments where only small Rm can be obtained. 
On the other hand, the gradual transition to MRI 
means that some of the typical features of the instabil- 
ity can be observed below the expected onset of the MRI. 

-Different structures for the MRI are obtained depend- 
ing on the hydrodynamical configuration. This is directly 
related to the imperfect nature of the bifurcation: the 
residual flow due to the boundaries strongly affects the 
geometry of the unstable mode. For instance, the mode 



6 



shown in Figure [3] has been obtained from the hydro- 
dynamical state given by eq. (|6]) and shown in Figure [TJ 
Without magnetic field, a weak outflow near the endcaps 
is produced, and the MRI mode generated in this case 
consists of poloidal vortices with the same outflowing jet 
near the endcaps. In other words, the hydrodynamical 
flow favors the upper branch of the pitchfork bifurcation, 
as in Figure HJ-bottom. A different hydrodynamical con- 
figuration would lead to a different MRI mode. Consider 
for instance the following rotation rates: 



Oi = 1. , tto = $h = ^4 = 0.13 



(10) 



Here, the endcap rings are rigidly rotating with the 
outer cylinder. The unmagnetized flow is characterized 
by a strong inflow near endcaps, and the outflow is now 
localized in the midplane. In the presence of a magnetic 
field, MRI modes show a similar geometry, with a strong 
outflow in the midplane. This corresponds to change 
h to —h in the equation (|9j), thus selecting the lower 
branch mode —A instead of A (this MRI mode is shown 
in the lower inset of Figure 2]) . 

-Finally, Figure HJ-bottom illustrates that in the case 
of an imperfect bifurcation, one of the branches is 
replaced by a saddle-node bifurcation of two modes, one 
stable and the other unstable, predicting the possible 
co-existence of two MRI modes with different geometries. 
Note that the transition to this new mode would be 
strongly subcritical. This is the equivalent of the well 
known " anomalous modes" observed in hydrodynamical 
Taylor- Couette flows. Due to its subcritical nature, this 
mode has not been observed in our simulations, the 
lower inset in Figure 2]-bottom corresponding to a sim- 
ulation with rotation rates given by eq. (fT0]h when the 
lower branch is continously connected to the Ekman flow. 



C. Saturation of the MRI 
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FIG. 5: Effect of the hydrodynamical Reynolds number 
Re on the growth of MRI modes. The magnetic 
Reynolds number and the Lundquist number are kept 
fixed to Rm = 30 and S = 4.5. Top: Time evolution of 

the radial magnetic field density A as a function of 
time, for different values of Re. Linear growth rates are 
independent of Re but the saturation level strongly 
depends on Re. Bottom: evolution of the saturated 
value of A as a function of Re. The saturation value of 
the MRI strongly decreases with Re. Red line indicates 
a fit corresponding to the scaling A ~ Re~ x l^ 



In all the simulations reported here, most of the pa- 
rameters have values comparable to those of laboratory 
experiments. For instance, in the Princeton experiment, 
the magnetic Reynolds number Rm can be varied from 
to 15, whereas the Lundquist number S is between and 
3. The exception is the hydrodynamical Reynolds num- 
ber, which can exceed 10 6 in the experiments, whereas 
our simulations are limited by numerical resolution to 
Re < 2 x 10 4 . It is thus important to understand how 
our numerical results scale to a system with Re ~ 10 6 . 
Figure [5] shows the evolution of the saturation of the MRI 
as a function of the Reynolds number. 

Figure [SJ-top shows the time evolution of the radial 
magnetic field density for different Reynolds numbers, 
ranging from Re = 10 3 to Re = 2 x 10 4 . In the linear 
phase, all the growth rates are the same, illustrating 
the relevance of the linear inviscid theory and the weak 
role played by the Reynolds number. On the contrary, 



the saturation level depends on the Reynolds number: 
as Re increases, the instability saturates to smaller 
values. It seems that the saturation of the MRI in our 
simulations scales like A ~ Re~^ , with f3 = 1/4. This 
is a dramatic illustration of the challenges of observing 
MRI experimentally, where Re must be several millions 
in order to achieve Rm > 1 because of the properties 
of liquid metals. However, a naive extrapolation of our 
numerical results suggests that the magnetorotational 
instability in the Princeton experiment should saturate 
at a value only one order of magnitude smaller than in 
the simulations reported here, and should be measurable. 

This scaling with Re is directly related to the nature of 
the non-linear saturation mechanisms involved in Taylor- 
Couette flows. Indeed, it has been suggested by Knobloch 
et al [20] that in such Taylor- Couette devices, saturation 
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mainly occurs by modification of the mean azimuthal pro- 
file: the shear in the flow is reduced, and is ultimately 
controlled by viscous coupling to the boundaries. Note 
that a simple dimensionless analysis shows that the ex- 
ponent j3 of our scaling corresponds to a balance between 
the Ekman pumping (~ Re -1 / 2 ) and the Maxwell stress 
tensor (~ B 2 ), in agreement with this scenario. However, 
it is important to note that the saturation of the MRI is 
expected to be strongly different in discs, where shear is 
fixed and boundaries may be stress- free. 



IV. NON-AXISYMMETRIC INSTABILITIES 

Let us now consider strongly magnetized situations. 
Figure [6] shows the same as Figure [TJ but for an applied 
field such that S = 7. In all the simulations reported 
here, the magnetic field is applied at the beginning of 
the simulation when the flow is at rest, letting the flow 
dynamically adapts as time evolves. At these large 
fields, we first note that the angular velocity profile is 
significantly different from the ideal Couette profile. 
In particular, the profile is steeper close to the inner 
cylinder and around r = (r\ + r2)/2. The distribution 
of the shear parameter q (Fig. [6fc) shows that the 
Stewartson shear layer (characterized by q > —3.5), 
confined to the endcaps in the purely hydrodynamical 
case, now penetrates deep inside the bulk of the flow 
[l7j : the applied magnetic field, by uniformizing the 
flow in the z-direction, reinforces and extends the shear 
between rings. This leads to the generation of a new 
free shear layer, different from the Stewartson layer. 
In addition, note that the poloidal flow is strongly 
affected by the applied field. The eight poloidal cells 
generated in the non-magnetic case are reinforced, and 
aligned in the z-direction. A strong axial jet develops at 
r — ( r i + r 2)/2, where the flow is ejected from the ring 
separation to the midplane. 

A similar magnetized free shear layer has been exten- 
sively described in spherical geometry. The shear layer 
appears when a strong magnetic field is applied to a 
spherical Couette flow, the magnetic tension coupling 
fluid elements together along the direction of magnetic 
field lines. In the case of an axial magnetic field, this cre- 
ates a particular surface E located on the tangent cylin- 
der, the cylindrical surface tangent to the inner sphere. 
Fluid inside E, coupled by the magnetic field only to the 
inner sphere, co-rotates with it, whereas fluid outside E 
(coupled to both spheres) rotates at an intermediate ve- 
locity. The jump of velocity on the surface E therefore 
results in a free shear layer, sometimes called Shercliff 
layer [21]. The thickness 7 of this MHD free shear layer 
is expected to vary like 7 ~ r/\/M, where M is the Hart- 
mann number. In our case, the Shercliff layer occurs on 
the cylindrical surface of radius r = (r\ + 7*2) /2, which 
separates regions of fluid coupled to the inner rings from 
fluid coupled to the outer rings. 




FIG. 6: Taylor- Couette flow when a strong uniform 
magnetic field B z is applied in the axial direction z, 

obtained for (^1,^3,^4,^2) = (400,180,65,53), 
Rm = 15, S = 7. Top-left: profile of Q at z = h/4. 
Top-right: Angular velocity £1 in (r, z) plane. We also 
show the profile obtained for S = 20. Bottom- left: U z in 

the (r, z) plane (streamlines indicate meridional flow. 
Bottom-right: Shear q in the (r, z) plane. The magnetic 

field tends to homogenize the flow in the z-direction. 
The jump of velocity between rings at endcaps extends 
into the bulk of the flow, generating a free shear layer at 
f — (^1 + r 2)/2. The 8 poloidal cells are re-enforced and 
a vertical jet is created. Note that this corresponds to 
the flow before the saturation of the non-axisymmetric 
instability. 



In spherical geometry, it has been shown that this flow 
configuration is unstable, and either the Shercliff layer or 
its associated poloidal flow leads to the generation of non- 
axisymmetric modes. In a recent work [9], it has been 
suggested that these modes, rather than MRI, were re- 
sponsible for the non-axisymmetric oscillations observed 
in the Maryland experiment. Interestingly, similar non- 
axisymmetric instabilities are observed in our cylindri- 
cal configuration, despite the difference in the geometry 
and the forcing. The red curve in Figure [2] is the locus 
of marginal stability for these non-axisymmetric modes; 
modes are unstable to the right (larger S for given Rm) of 
the curve. The critical Reynolds number seems to follow 
the scaling Rm c ~ S 2 . These non-axisymmetric modes 
are also restabilized for strong values of the applied field, 
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Time [s] 

FIG. 7: Time evolution of Fourier components of the 
azimuthal velocity. The shear layer is destabilized to 
non-axisymmetric modes via a Kelvin-Helmoltz type 
instability. In the linear phase, the m = 3 azimuthal 
harmonic is the most unstable; the saturated state is 
dominated by m = 1. Inset: the same evolution as a 
semilog plot. 



but we do not have enough points to determine the cor- 
responding scaling. 

The marginal-stability curve Rm c ~ S 2 corresponds 
in fact to the region of the parameter space where the 
Elsasser number, A = B 2 /(fjbprjQ), is equal to unity. 
From Figure [2j it can be seen that our shear layer 
instability extends to small magnetic Reynolds number 
(Rm < 1). Apparently, the shear layer instability 
reported here is induct ionless, meaning that the time 
derivative of the magnetic field can be neglected in 
the induction equation ([3]). This is a crucial difference 
from the standard magnetorotational instability, for 
which induction is essential (Rm > 1). The fluid 
Reynolds number, Re, is more relevant to the present 
instability than Rm. This induct ionless property is also 
possessed by the so-called Helical MRI, which has been 
understood as an inertial wave weakly destabilized by 
the magnetic field [22|, [23[. The inductionless nature 
of the shear-layer instability is easy to understand 
physically: although a strong magnetic field is necessary 
to generate the layer, the destabilization of this layer is 
triggered by a hydrodynamical shear instability, similar 
to Kelvin-Helmoltz, whose form is nearly constant along 
the field lines and therefore scarcely affected by magnetic 
tension. The instability is suppressed if the difference 
between the rotation rate of the two rings is too small. 

Figure [71 shows the time evolution of azimuthal Fourier 
modes for Rm = 15 and S = 11. The initially 
axisymmetric MHD shear layer is destabilized to sev- 
eral azimuthal wave numbers. In the linear phase, the 



m = 3 mode clearly dominates. As the instability sat- 
urates, there is a cascade of energy towards lower az- 
imuthal wavenumbers, and the nonlinear saturated state 
is strongly dominated by an oscillating m = 1 mode. 
In the saturation phase, the thickness of the shear layer 
is increased, indicating that the instability saturates by 
modifying the shear profile. 

Figure [8] shows the structure of these non-axisymmetric 
instabilities. In the linear phase, the structure is very 
similar to Kelvin-Helmoltz modes of the free shear layer, 
taking the form of columnar vortices transverse to the 
(r, (/)) plane, independent of the axial direction and 
therefore symmetrical with respect to the equatorial 
plane. Note that the vortices spiral slightly in the 
prograde direction, as expected from a Kelvin-Helmoltz 
instability in cylindrical geometry transporting the 
angular momentum outward. Figure [SJ-top shows an 
isosurface of 25% of the total shear q in the flow. This 
illustrates how the rotational symmetry of the free 
shear layer is broken in the <\> direction, the initially 
axisymmetric layer rolling up in a series of horizontal 
vortices by a mechanism similar to Kelvin-Helmoltz 
instabilities. 

The structure of the magnetized shear layer strongly 
depends on the parameters. For a fixed value of the 
global rotation, the thickness I of this shear layer scales 
with the Hartmann number as I ~ r/y/~M, a scaling sim- 
ilar to what is observed in spherical geometry. On the 
other hand, if the applied field is decreased after this 
shear layer has formed, the vertical extent of the shear 
layer is reduced, and eventually reaches a critical length 
where the Kelvin-Helmoltz instability does not occur. In 
our simulations, it appears that the vertical extent of the 
shear layer scales as h ~ rM/y/Re = ry/X. The condi- 
tion h ~ r thus gives the A = 1, destabilization condi- 
tion observed in our simulations: the non-axisymmetric 
modes develop by Kelvin-Helmoltz destabilization of the 
shear layer, provided that the layer is sufficiently extended 
in the vertical direction (A > 1). 



V. COMPARISON WITH LABORATORY 
EXPERIMENTS 

Finally, we compare our numerical simulations to re- 
sults obtained in the Princeton MRI experiment, where 
a similar setup is used (see [24[ for more details on the 
experimental results mentioned here). Because of the 
very small magnetic Prandtl number of liquid metals 
(Pm — Rm/Re < 10 -5 ), it is very challenging to achieve 
fluid Reynolds numbers large enough to observe MRI in 
liquid- metal experiments. In Figure [2j the region of pa- 
rameter space accessible to the Princeton experiment is 
indicated by the dashed square. 

Our numerical simulations indicate that the magne- 
torotational instability should in theory be observable in 
the Princeton MRI experiment between Rm = 10 and 
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FIG. 8: Structure of the non-axisymmetric instability 
for Rm = 15, S = 11. Top: Structure in the r, </> plane, 
for z = H/A. Non-axisymmetric component of the 
radial velocity u r in the equatorial plane during the 
linear phase. Bottom: 3D isosurface of the shear q. The 
mode consists of non-axisymmetric axial columnar 
modes, i.e. vortices in the r, (f) plane (the axial 
component is much weaker) and mainly independent of 
the z-direction. Isosurface of the shear shows that these 
modes are essentially Kelvin-Helmoltz destabilization of 
the central magnetized shear layer. 



Rm = 15; the latter corresponds to the maximum de- 
signed rotation rotation speed. Up to now, the maximum 
magnetic Reynolds number reached in the Princeton ex- 
periment is Rm = 9, i.e. just below threshold. Our 
simulations suggest that the MRI should take the form 
of a gradual amplification of the residual Ekman flow. 
Since the experiment is still operating very close to the 
expected onset of the MRI, it is necessary to operate at 
higher Rm to see MRI clearly. A further difficulty is that 
the saturation level is expected to be smaller in the ex- 
periment that in the simulations because of the Re~ 1 ^ 
scaling discussed above (Fig. [5]). 

By contrast, the shear layer instability is inductionless 
and should be observed at much smaller Rm. Figure [9] 
compares the numerical and experimental results for 
the horizontal structure of the mode. Non-axisymmetric 
modes have indeed been observed at very low rotation 
rates in the experiment, and these modes present several 
features identical to the instabilities reported here: 



• First, the generation of non-axisymmetric modes 
in the experiment follows the marginal stability curve 
A c = 1, in agreement with our numerical prediction. 
Therefore, the experimental instabilities are also in- 
ductionless: they have been observed for rotation rates 
of the inner cylinder as small as Qi = 0.05 rads -1 
(Rm <C 1). 

• The structure of the modes is very similar. Figure [9] 
compares the geometry of the modes in the (r, (/)) plane 
obtained in the experiment (bottom) and our simulations 
(top), for a case in which the dimensionless experimental 
parameters match the numerical ones except in mag- 
netic Reynolds number. The similarity of the instability 
across two orders of magnitude in Rm confirms its 
inductionless nature. In both cases, non-axisymmetric 
modes consist of equatorially symmetric vortices in the 
(r, z) plane taking the form of an m = 1 spiral mode 
sheared in the azimuthal direction. Note however that 
in the simulation, close to the instability onset, the 
modes spiral retrogradely in the midplane, but tend to 
be more prograde close to the rings (in agreement with 
the expected structure of a Kelvin-Helmoltz instability). 
This effect has not been observed in the experiment, 
where the spiral seems to be retrograde even relatively 
close to the endcaps. This could be due to the effect 
of the mean shear flow, or the secondary excitation of 
typical waves of the system, such as magneto-coriolis 
waves [6[. A prograde spiral is however expected 
sufficiently close to the endcaps, since it corresponds to 
an outward transport of the angular momentum by the 
Kelvin-Helmoltz instability. 

• The time evolution is similar. In both numerical 
simulations and experiment, the saturated state is 
an m = 1 mode rotating at a constant speed in the 
equatorial plane. Before reaching this m = 1 mode, the 
experimental flow exhibits transient oscillations with 
different azimuthal symmetry, namely an m = 3 followed 
by an m = 2 mode. This is similar to the behavior of 
our numerical model, reported in Figure [71 

• In the experiment, non-axisymmetric modes are ob- 
served only when the difference of angular velocity be- 
tween the inner ring and the outer ring is sufficiently 
large, and global rotation can stabilize the flow. These 
two features are also observed in our numerical simula- 
tions. 

These several similarities strongly suggest that the 
non-axisymmetric modes observed in the Princeton MRI 
experiment are of the same nature as the ones reported 
in our numerical simulations, namely Kelvin-Helmoltz- 
type destabilization of a free shear layer created by the 
conjugate action of the axial magnetic field and the 
jump of angular velocity between the rotating rings at 
the endcaps. Note however that the magnetic boundary 
conditions at the endcaps are not exactly the same: our 
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FIG. 9: Comparison between numerical simulation 
(top) and ultrasound measurements Q from the 
Princeton experiment (bottom). In both cases, inner 
rings corotate with the inner cylinder (Qi = ^3), outer 
rings and cylinder are at rest (Q2 = = ^4), 
Re = 1450, and A = 1.6. Simulation has Rm = 0.02, 
and experiment, Rm = 0.002. Figures show the 
non-axisymmetric component of the azimuthal velocity 
in the midplane, when the Shercliff layer is unstable. 
A very good agreement is obtained: in both cases, the 

saturated state consists of an m = 1 spiral vortices 
rotating in the azimuthal direction. In the experiment, 
Uff, is obtained from two ultrasound probes measuring 
the radial profile of the velocity at two different 
azimuthal positions. 



simulations use the idealized 'pseudo-vaccum' boundary 
conditions given by eq.flSJ), whereas the endwalls are 
insulating in the Princeton experiment. Since insulating 
endwalls do not couple with magnetic field, the role 
of the electrical currents close to the endcaps in the 
establishment of the layer could be slightly different 
between experiments and numerics. 

Interestingly, similar instabilities occur in spherical ge- 
ometry, and could be involved in the oscillations observed 
in the Maryland experiment. More recently, similar ro- 
tating rings have also been used in the PROMISE experi- 
ment in order to obtain clearer observation of the Helical 
MRI, and non-axisymmetric m = 1 modes have been 
reported. It would be interesting to see how these oscil- 
lations compare to the non-axisymmetric modes reported 



here. 



VI. CONCLUSION 

In this article, we have reported three-dimensional nu- 
merical simulations of a magnetized cylindrical Taylor- 
Couette flow inspired by laboratory experiment aiming 
to observe the magnetorotational instability. 

The influence of the boundaries on the MRI has been 
studied. We have shown that that the finite axial extent 
of the flow deeply modifies the nature of the bifurcation. 
In the presence of realistic boundaries, the MRI appears 
from a strongly imperfect supercritical pitchfork bifurca- 
tion, leading to several predictions for the observation of 
the instability in the laboratory. First, the transition to 
MRI is expected to gradually emerge from the residual 
recirculation driven by the top and bottom boundaries. 
Moreover, this Ekman flow directly modifies the struc- 
ture of the instability by constraining the geometry of 
the mode, and different MRI modes can be observed de- 
pending on the hydrodynamical background state. As a 
consequence, the endcap rings of the Princeton experi- 
ment are not only useful to reduce the Ekman recircula- 
tion, but they can also be used to select different MRI 
modes. 

In the second part of this study, we have reported the 
generation of non-axisymmetric modes when the applied 
magnetic field is sufficiently strong compared to the rota- 
tion (Elsasser number > 1). Although they present some 
similarities with the MRI, these modes are of a very dif- 
ferent nature: when the axial magnetic field is sufficiently 
strong, it generates a new free shear layer in the flow, by 
extending the jump of angular velocity between the inner 
and outer endcap rings into the bulk of the flow. When 
the shear is sufficiently strong, this layer is destabilized 
by non-axisymmetric modes of a Kelvin-Helmoltz type. 
It is interesting to note that this new instability shares 
several properties with the MRI: It transports angular 
momentum outward, it is stabilized at strong applied 
field, and a critical magnetic field is needed. However, 
these modes are inductionless, and therefore extend to 
very small Rm. This is an important difference from the 
classical MRI. 

To finish, we have compared our numerical simulations 
to experimental results from the Princeton experiment. 
Good agreement is obtained at small Rm, where non- 
axisymmetric oscillations are also observed in the exper- 
iment. The marginal stability curve, the geometry of the 
mode, and the time evolution of these non-axisymmetric 
oscillations are similar to our numerical simulations. This 
indicates that the non-axisymmetric oscillations observed 
in the Princeton experiment could be related to the desta- 
bilization of the free shear layer (Shercliff layer) reported 
in the present numerical work. 

To summarize, the present numerical study shows that 
the generation, structure, and saturation of the MRI are 
strongly modified by the presence of no-slip boundaries. 
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In addition, new instabilities, similar to the MRI but 
inductionless, are generated in the presence of a strong 
magnetic field. This suggests that the observation of the 
magnetorotational instability in a laboratory experiment 
could be radically different from what is expected from 
local theory, or even axially-periodic Taylor- Couette sim- 
ulations. However, careful comparisons between numer- 
ical simulations and experimental results could lead to 
a clear identification of the magnetorotational instability 
in the laboratory, and to a better understanding of the 
angular momentum transport. 
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